IMPLICIT INTEGRATION OF THE TIME-DEPENDENT 
GINZBURG LANDAU EQUATIONS OF SUPERCONDUCTIVITY 

D. O. GUNTER , H. G. KAPER , AND G. K. LEAF* 

Abstract. This article is concerned with the integration of the time-dependent Ginzburg- 
Landau (TDGL) equations of superconductivity. Four algorithms, ranging from fully explicit to fully 
implicit, are presented and evaluated for stability, accuracy, and compute time. The benchmark 
problem for the evaluation is the equilibration of a vortex configuration in a superconductor that is 
embedded in a thin insulator and subject to an applied magnetic field. 
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1. Introduction. At the macroscopic level, the state of a superconductor can 
be described in terms of a complex- valued order parameter and a real vector potential. 
These variables, which determine the superconducting and electromagnetic properties 
of the system at equilibrium, are found as solutions of the Ginzburg-Landau (GL) 
equations of superconductivity. They correspond to critical points of the GL energy 
functional Jy, [2J, so in principle they can be determined by minimizing a functional. 
In practice, one introduces a time-like variable and computes equilibrium states by 
integrating the time-dependent Ginzburg-Landau (TDGL) equations. The TDGL 
equations, first formulated by Schmid |3| and subsequently derived from microscopic 
principles by Gor'kov and Eliashberg |ij, are nontrivial generalizations of the (time- 
independent) GL equations, because the time rate of change must be introduced in 
such a manner that gauge invariance is preserved at all times. 

\n ■ We are interested, in particular, in vortex solutions of the GL equations. These 

are singular solutions, where the phase of the order parameter changes by 2ir along 

0"N ■ any closed contour surrounding a vortex point. Vortices are of critical importance in 

technological applications of superconductivity. 
f~| ■ Computing vortex solutions of the GL equations by integrating the TDGL equa- 

tions to equilibrium has the advantage that the solutions thus found are stable. At 
the same time, one obtains information about the transient behavior of the system. 
Integrating the TDGL equations to equilibrium is, however, a time-consuming pro- 
cess requiring considerable computing resources. In simulations of vortex dynamics in 
superconductors, which were performed on an IBM SP with tens of processors in par- 
allel, using a simple one-step Euler integration procedure, we routinely experienced 
equilibration times on the order of one hundred hours H, R, ■ Incremental changes 
would gradually drive the system to lower energy levels. These very long equilibration 
times arise, of course, because we are dealing with large physical systems undergoing a 
phase transition. The energy landscape for such systems is a broad, gently undulating 
plain with many shallow local minima. It is therefore important to develop efficient 
integration techniques that remain stable and accurate as the time step increases. 

In this article we present four integration techniques ranging from fully explicit 
to fully implicit for problems on rectangular domains in two dimensions. These two- 
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dimensional domains should be viewed as cross sections of three-dimensional systems 
that are infinite and homogeneous in the third direction (orthogonal to the plane of 
the cross section), which is the direction of the field. The algorithms are scalable 
in a multiprocessing environment and generalize to three dimensions. We evaluate 
the performance of each algorithm on the same benchmark problem, namely, the 
equilibration of a vortex configuration in a system consisting of a superconducting 
core embedded in a blanket of insulating material (air) and undergoing a transition 
from the Meissner state to the vortex state under the influence of an externally applied 
magnetic field. We determine the maximum allowable time step for stability, the 
number of time steps needed to reach the equilibrium configuration, and the CPU 
cost per time step. 

Different algorithms correspond to different dynamics through state space, so the 
eventual equilibrium vortex configuration may differ from one algorithm to another. 
Hence, once we have the equilibrium configurations, we need some measure to assess 
their accuracy. For this purpose we use three parameters: the number of vortices, 
the mean intervortex distance (bond length), and the mean bond angle taken over 
nearest-neighbor pairs of bonds. When each of these parameters differs less than a 
specified tolerance, we say that the corresponding vortex configurations are the same. 

Our investigations show that one can increase the time step by almost two orders 
of magnitude, without losing stability, by going from the fully explicit to the fully 
implicit algorithm. The fully implicit algorithm has a higher cost per time step, but 
the wall clock time needed to compute the equilibrium solution (the most important 
measure for practical purposes) is still significantly less. All algorithms yield the same 
equilibrium vortex configuration. 

In Section |2j, we present the Ginzburg-Landau model of superconductivity, first 
in its formulation as a system of partial differential equations, then as a system of 
ordinary differential equations after the spatial variations have been approximated by 
finite differences. In Section y, we give four algorithms to integrate the system of 
ordinary equations: a fully explicit, a semi-implit, an implicit, and a fully implicit 
algorithm. In Section [|, we present and evaluate the results of the investigation. The 
conclusions are summarized in Section pL 

2. Ginzburg-Landau Model. The time-dependent Ginzburg-Landau (TDGL) 
equations of superconductivity |2j, y, ^ are two coupled partial differential equations 
for the complex- valued order parameter i/j = l^e 1 ^ and the real vector- valued vector 
potential A, 

<"> 2^0 (I + T # ) * -i (? v - 7 A ) 2 * + "* - «'' a *' 

(2.2) „(___+v*]= VxVxA + J„ 

\c at J Att 

Here, J s is the supercurrent density, which is a nonlinear function of "0 and A, 

(2.3) J 8 = ^(0*W ~ W*) - ^HVfA = fL^f (hVcb - ^A) . 

lim s m s c m s \ c / 

The real scalar-valued electric potential $ is a diagnostic variable. The constants in 
the equations are K, Planck's constant divided by 2tt; a and b, two positive constants; 
c, the speed of light; m s and e s , the effective mass and charge, respectively, of the 
superconducting charge carriers (Cooper pairs); v, the electrical conductivity; and 
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D, the diffusion coefficient. As usual, i is the imaginary unit, and * denotes complex 
conjugation. 

The quantity |i/'| 2 represents the local density of Cooper pairs. The local time 
rate of change dt A of A determines the electric field, E = (f /c)dt A + V<f>, its spatial 
variation the (induced) magnetic field, B = V x A. 

The TDGL equations describe the gradient flow for the Ginzburg-Landau energy, 
which is the sum of the kinetic energy, the condensation energy, and the field energy, 



(2.4) E = 



I 



2m s 



A V 



+ -«M + dff + |VxA| 



dx. 



A thermodynamic equilibrium configuration corresponds to a minimum of E. 



The energy functional (2.4) assumes that there are no defects in the supercon- 
ductor. Material defects can be naturally present or artifically induced and can be in 
the form of point, planar, or columnar defects (quenched disorder). A material defect 
results in a local reduction of the depth of the well of the condensation energy. A 
simple way to include material defects in the Ginzburg-Landau model is by assuming 
that the parameter a depends on position and has a smaller value wherever a defect 
is present. 

2.1. Dimensionless Form. Let tp 2 ^ = a/b, and let A, £, and H c denote the 
London penetration depth, the coherence length, and the thermodynamic critical 
field, respectively, 



(2.5) 



/ 2 \ 1/2 

/ m s c \ 



X {Wl>e* ' 



/ p, 2 \ 1/2 

(£-) , H c = ( 47 ra^) 1/2 - 
\2m s a J 



In this study, we render the TDGL equations dimensionless by measuring lengths in 
units of £, time in units of the relaxation time £ 2 / ' D, fields in units of H c ^/2, and 
energy densities in units of (l/Air)H 2 . The nondimensional TDGL equations are 



(2.6) 

(2.7) 



— +i$W=(V-~A) V + rip - IV'lV) 

dA _ _ _ 

a \ — + kV$ ) = -V x V x A + J s 



where 



(2. 






V0--A 

K 



Here, k = \/£, is the Ginzburg-Landau parameter and a is a dimensionless resistivity, 
a = (A-kD/c 2 )v. The coefficient r has been inserted to account for defects; r(x) < 1 if 
x is in a defective region; otherwise t(x) = 1. The nondimensional TDGL equations 
are associated with the dimensionless energy functional 



(2.9) 



E 



V 



A V 



+ (-r|^| 2 + i|0| 4 ) + |VxA| 2 



da;. 
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2.2. Gauge Choice. The (nondimcnsional) TDGL equations are invariant un- 
der a gauge transformation, 



(2.10) 



Q x : (V, A, $) .- tye* A + kV X , $ ~ d tX ). 



Here, x can be any real scalar-valued function of position and time. We maintain the 
zero-electric potential gauge, $ = 0, at all times, using the link variable U, 



(2.11) 



U = exp 



- / A 

K 



This definition is componentwise: U x = cxp(— in J A x (x' ,y, z) dx'), 
gauged TDGL equations can now be written in the form 



The 



(2.12) 


fi—x,y.z 


(2.13) 


dA. 

a— =-VxVxA + J s , 
at 


where 




(2.14) 


Js,a = - Im 
K 


(u^r^u^) 


, n = x,y,z 



2.3. Two-Dimensional Problems. From here on we restrict the discussion to 
problems on a two-dimensional rectangular domain (coordinates x and y), assuming 
boundedness in the x direction and periodicity in the y direction. The domain rep- 
resents a superconducting core surrounded by a blanket of insulating material (air) 
or a normal metal. The order parameter vanishes outside the superconductor, and 
no superconducting charge carriers leave the superconductor. The whole system is 
driven by a time-independent externally applied magnetic field H that is parallel to 
the z axis, H = (0, 0, H). The vector potential and the supercurrent have two nonzero 
components, A = (A Xl A y ,Q) and J s = (J x , J y ,0), while the magnetic field has only 
one nonzero component, B = (0, 0, B), where B = d x A y — d y A x . 

2.4. Spatial Discretization. The physical configuration to be modeled (su- 
perconductor embedded in blanket material) is periodic in y and bounded in x. In 
the x direction, we distinguish three subdomains: an interior subdomain occupied 
by the superconducting material and two subdomains, one on either side, occupied 
by the blanket material. We take the two blanket layers to be equally thick, but do 
not assume that the problem is symmetric around the midplane. (Possible sources 
of asymmetry are material defects in the system, surface currents, and different field 
strengths on the two outer surfaces.) 

We impose a regular grid with mesh widths h x and h y , 

(2.15) ttij = {xuXi+i) x (yj,y j+ i), x { = x + ih x ; y = y Q + jh y , 



assuming the following correspondences: 

Left outer surface: x = Xo 
Left interface: x = x ns 
Right interface: x = x ric 
Right outer surface: x = x Uq 



2 ,ix 



2 x 

i/l 

2 IL x-) 



2 x 



o. 
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One period in the y direction is covered by the points j — 1, . . . ,n y . We use the 
symbols Sc and Bl to denote the index sets for the superconducting and blanket 
region, respectively, 

(2.16) Sc = {{i,j) : (i,j) e [n sx ,n ex ] x [l,n y ]}, 

(2.17) Bl = {(i,j) : (i,j)£ [l,n ax -l]U[n e x + l > n B ] x [1,%]}. 

The order parameter ip is evaluated at the grid vertices, 

(2.18) ip itj = f(xi,yj), (i,j) E Sc, 

the components A x and A y of the vector potential at the midpoints of the respective 
edges, 

(2.19) A x - itj =A x (xi + ±h x ,yj), A y . itj = A y {x i ,y j + ^h y ), (i,j)eScUBl, 
and the induced magnetic field B at the center of a grid cell, 



(2.20) 



B, 



B(x. t + \h Xl yj + \h y ) 

-<*y;i-{-l,j ~ -<*-y;i.j s±xyi,j-\-l 



A 



XVI,] 



^ (»,j)6ScUBl, 



see Fig. 2.1. The values of the link variables and the supercurrent are computed from 

J+1 4 x i 



A y;U X 




Fig. 2.1. Computational cell with evaluation points for ip, A x , and A y . 

the expressions 

('2 211 U ■ ■ — p- iK ~ lh * A *;i,J TJ . . — p -iK~ 1 hyAy,i }j 

(2.22) J x;it j = -r-Im [ipijUxfrjipi+ij], Jy,i,j = — j— Im [iptjUyjjipij+i] 



kHz 



The discrctizcd TDGL equations are 



t\iflq. 



(2.23) 

(2.24) a- 

(2.25) a 



dt 



dipjj 

h 

dA x -i 

~~dT 

dt 



(L xa (U Xi .^)il).j) i + (Ly V (Uy.i,.)il>i,.) i + N(il)i d ) i (i,j) e Sc, 



\ W x\i,-)j \*-'yx-t^ 



yx-™-v,-,-Ji,j 



Jx;i,j, (i,j)GScUBl, 



{D xx A y .., j ) i - (D xy A x ..,.) +J y . tJ . (i,j) e ScUBl, 
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where 



(2.26) {L xx (U x .. d )ip.j) t = h~ 2 [U^ijipi+ij - 2V»ij + 17*. i _ 1)j i/'«-i,j] , 

(2.27) {L V y{U v . it .)ip ii .) j = hy 2 [Uy-ijtpij+x - 2^ + t^-iV'i.j-i] , 
(2-28) TV (Viy) = nj-ipij - \ipij |Vy , 

(2.29) (-DwAr-A.)j = V 2 [4rAi+l ~ 2A *;« + ArAj-l] , 

(2.30) (D XX Ay.., j ) i = ft" [A.y ; j + lj - 2A y;ij + Ay ;i _lj] , 

(2.31) (PvxAy,-,-) it j=h~ h y [(Ay ii+li j — Ay-ij) — (Ay-i + lj-X — Ayfrj-l) 

(2.32) \DxyAx;-,-)ij — h x h y [(A x -ij + i — ^x;i,j) ~ {A x ;i-l,j+l — Ax;i-\,j , 

The interface conditions are 

(2.33) i> nsx -l,j = Ux-n^-ijlpn^j, i„+ij = U*;n ex ,ji J n ex ,j, j = 1, . . . ,n y . 

At the outer boundary, B is given, 

(2.34) Bo,, = H Lj , B„^ =H Ri ,j=l,..., n v . 
The resulting approximation is second-order accurate . 



3. Time Integration. We now address the integration of Eqs. (2.23)— (p.25) 



The first equation, which controls the evolution of t[>, involves the second-order linear 
finite-difference operators L xx and L yy , whose coefficients depend on A x and A y , 
and the local nonlinear operator N, which involves neither A x nor A y . Each of the 
other two equations, which control the evolution of A x and A y respectively, involves 
likewise a second-order linear finite-difference operator, but with constant coefficients, 
and the nonlinear supercurrent operator, which involves ip, A x , and A y . The following 
algorithms are distinguished by whether the various operators are treated explicitly 
or implicitly. 

3.1. Fully Explicit Integration. Algorithm I uses a fully explicit forward 
Euler time-marching procedure for ip, A x , and A y . Starting from an initial triple 
(ip ,A®,A°), we solve for n = 0,1,..., 

(3 1} tiA a ~ Ei = {L aa (JJ2.. ii Wj) i + (L yy {ui it Wi). + N (^) , (i t j) € Sc, 

A n+1 — A n 
(3.2) a *™ *Jhi = (D vy Al. it ) - (D yx A n y ,. t ). + J x " ;lJ , (i,j) e ScU Bl, 



At 

An+l _ An 
A,.," , A y-i 



( 3.3) a v™ yjid = (D xx A n y ,. t] ). - (D xy A2.. t .). + J^. (i,j) e ScU Bl, 



where J" is defined in terms of ip n , A™ , and A y in the obvious way. The initial triple 
is usually chosen so the superconductor is in the Meissner state, with a seed present 
to trigger the transition to the vortex state. 

Algorithm I has been described in M. It has been implemented in a distributed- 
memory multiprocessor environment (IBM SP2); the transformations necessary to 
achieve the parallelism have been described in || . The code uses the Message Passing 
Interface (MPI) standard [[LO] as implemented in the MPICH software library |ll[] 
for domain decomposition, interprocessor communication, and file I/O. The code has 
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been used extensively to study vortex dynamics in superconducting media || f| [/J. 
The underlying algorithm provides highly accurate solutions but requires a significant 
number of time steps for equilibration. For stability reasons, the time step At cannot 
exceed 0.0025. 

3.2. Semi-Implicit Integration. Algorithm II is generated by an implicit treat- 
ment of the second-order linear hnite-difference operators D yy and D xx in the equa- 
tions for A x and A y , respectively, 

(3 4) tiJ ^ i%i = {L wa (P5 i . J )1>?j) i + {LyyiUl^l). + N (V&) , (i,j) G Sc, 

A n+l — A n 

(3.5) a S1M_ *Jhl = {D yy Al£). - {D yx AU . . + J™^, (i,j) G Sc U Bl, 

A n+1 — A n 

( 3. 6) a ^L_y±L = (D XX A^). - (D xy A2.. t .). t . + J^ y (i,j) G ScU Bl. 

Equations fl3.5J ) and ( |3.6| ) lead to two linear systems of equations, 

(3.7) {l-^D w \A^ = Fi^ n ,A^,A^), i = 1,... ,n x - 1, 

1 ■• ' 1™+! — n.(„i, n A n A n \ 
y;j 



a 



(3.8) [I-—D xx )A^ = G j (r,A:,A^), j = l,...,n y , 



for the vectors of unknowns A x -i — {A x -ij : j = 1,. . . ,n y } and j4j, ; j = {A y ^_j : i = 
1, . . . ,n x — 1}. The matrix D yy has dimension n y x n y and is periodic tridiagonal 
with elements —h~ 2 , 2h~ 2 , —h~ 2 ; the matrix D xx has dimension (n x — 1) x (n x — 1) 
and is tridiagonal with elements ~h x 2 , 2h x 2 , —h x 2 , (except along the edges, because 
of the boundary conditions). Both matrices are independent of i and j. Furthermore, 
if the boundary conditions are time independent, they are cons tant throu gho ut the 



time-stepping process. Hence, the coefficient matrices in Eqs. (3/7) and (3J3) need 
to be factored only once; in fact, the factorization can be done in the preprocessing 
stage and the factors can be stored. 

In a parallel processing environment, the coefficient matrices extend over several 



processors, so Eqs. (3.7) and (3.£) are broken up in blocks corresponding to the manner 
in which the computational mesh is distributed among the processor set. We first solve 
the equations within each processor (inner iterations) and then couple the solutions 
across processor boundaries (outer iterations). Hence, we deal with interprocessor 
coupling in an iterative fashion. Two to three inner iterations usually suffice to reach 
a desired tolerance for convergence. After each inner iteration, each processor shares 
boundary data with its neighbors through MPI calls. 

3.3. Implicit Integration. Algorithm III combines the semi-implicit treatment 
of A x and A y with an implicit treatment of the order parameter, 

(3 9) ^ ^ ^ = (L ax (US.. d )^ 1 ) i + (L yy (U- l ^l +l ) j + N ty&) , (i,j) G Sc, 

A n+1 — A n 
(3.10)7 x ™ ^ X ' M = (DyyAZg). - (D yx A^.. t .). t . + J x " ;iji , (i, j) G Sc U Bl, 

A n+1 — A n 
(3.11V ^ At ™ = (D XX A^). - (D xv A2,,.). t . + J%. id . (i,j) G ScU Bl. 
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The second and third equation are solved as in the semi-implicit algorithm of the 
preceding section. The first equation is solved by a method similar to the method of 
Douglas and Gunn |L2] for the Laplacian. 



We begin by transforming Eq. (3.9) into an equation for the correction matrix 
t> n+l = ip n+1 — tb n . The equation has the general form 



(3.12) 



(/ - At(L xx + L yy )) <b n+l = F{r,A n x , A n y ) 



If At is sufficiently small, we may replace the operator in the left member by an 
approximate factorization, 



(3.13) (/ - At(L xx + L yy )) « (I - AtL xx ) (I - AtL yy ) , 

and consider, instead of Eq. (BT2J), 



(3.14) 



(I - AtL XX ) (I - AtLyy) 



/,«+! 



F(r,A n x ,A"). 



This equation can be solved in two steps, 



(3.15) 
(3.16) 



(/ - AtL xx ) <p = F, 
ip. 



(I - AtLyy) fi^ 1 



The conditions (2.33), which must be satisfied at the interface between the supercon- 
ductor and the blanket material, require some care. If we impose the conditions at 
every time step, then 



^riax—l,] ~ ^x;n sx -l,j'Pn ea .,j + Wx;n etc -l,j ^x;n ex -l,j\ V* 

(TI n+1 )* - (U n V 

V x;n e3 .,j J \ x;n sx — l,j / 



n B x:,3> 



in+1 _ (Tjn+1 )* An+1 

y n cx + lj — \ u x;n ex ,j) Vn ex ,j 



K 



for j = 1, . . . ,n y . These conditions couple the correction <j> to the update of A x . To 
eliminate this coupling, we solve Eq. (3.1 2|) subject to the reduced interface conditions 



(3.17) 
(3.18) 



b n+1 — U n+1 S n+1 i — 1 n 



+1 = (U n+l )* A n+1 7 = 1 n 



When Eq. (3.12) is replaced by Eq. (3.14), these conditions are inherited by the 
system (3.15). 



3.4. Fully Implicit Integration. Algorithm IV uses a fully implicit integration 
procedure for the order parameter, 

(3.19) ^ ^ ^ = {L xx (U2.. J W$ 1 ) i + (L VV (U^.)^; 1 ). + N (^ +1 ) , (i,j) e Sc, 

A n+l — A n 
(3.20)7 X ^ M x ^ = {DyyA n x + 1 ) j - (D vx A^,). tj + Jl iid , (i,j) G Sc U Bl, 

A n+1 — A n 
(3.21)0- mi ' j . , mi ' j = (D x , r A n +]).- (I), „,!'.:.. 



At 



y;-,j/i 



'xy^ x - 



Jy;i,j- (i,j)eScUBL 



The new element here is the term N (ipff ) in the first equation. 
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The second and third equations are solved again as in the semi-implicit algorithm. 
The first equation is solved by a slight modification of the method used in the im- 
plicit algorithm of the preceding section, The modification is brought about by the 
approximation 

(3.22) N (ip n+1 ) = Tijj n+1 - \ip n+1 \ 2 ip n+1 ~-^-(S (tp n ) - 4> n ) , 

where S is a nonlinear map, 

T 1 ^ 



(3.23) S{ip) 



[H 2 + (T-|Vf)exp(-2TAi)] 



1/2- 



(This approximation is explained in the remark below.) Equation (3.19) is again of 



the form (3.12), but with a different right-hand side, 

(3.24) (7 - At(L xx + L yy M l+1 - G(^\ A n xl A n y ). 



The difference is that, where F in Eq. (3.1 2|) contains a term (At)N (ip n ), G in 



Eq. (|3.24) contains the more complicated term S (ip n ) — ip n - 



Remark. The approximation (3.22) is suggested by semigroup theory. Symboli- 
cally, 

(3.25) N&) = lim 5(A y~^ . 

v ' v ' At^O At 

To find an expression for the "semigroup" S, we start from the continuous TDGL 



equations (J2.6J) (2.8 ) (zero-electric potential gauge, $ = 0), using the polar represen- 
tation tp = \ip\e^, 

(3.26) d t \il>\ = A|V| - M|V0 - n~ l A\ 2 + T \il>\ - M 3 , 

(3.27) \tp\d t (J3 = 2(V]VI) ' (V</> - kT x A) + \ip\V ■ (V0 - k _1 A), 

(3.28) ad t A = -V x V x A + k- x |^| 2 (V0 - k -1 A). 

At this point, we are interested in the effect of the nonlinear term |^| 3 on the dynamics. 
To highlight this effect, we concentrate on the time evolution of the scalar u = \i/)\ and 
the vector v = V</> — k _1 A. (In physical terms, u 2 is the density of superconducting 
charge carriers, while u 2 v is k times the supercurrent density.) Ignoring their spatial 
variations, we have a dynamical system, 

(3.29) u = -u\v\ 2 + tu - u 3 , 

(3.30) v' = -eu 2 v, 

where ' denotes differentiation with respect to t, and e = (k 2 ct) _1 . This system yields 
a pair of ordinary differential equations for the scalars x = u 2 and y = \v\ 2 , 

(3.31) x' = 2x{r-x-y), 

(3.32) y' = -2exy. 

If k is large, e is small, and the dynamics are readily analyzed. To leading order, y 
is constant; y = is the only meaningful choice. (Recall that xy 1 / 2 is k times the 
magnitude of the supercurrent density.) Then the dynamics of x are given by 

(3.33) x'=2x(t-x). 
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We integrate this equation from t = t n to t, 

(3-34) x(t) = —— s -. J ,, . - . r-. 

y ' x(t n ) + (t - x(t n )) exp(-2r(t - £„)) 

In particular, 

(3.35) x(t n+1 ) - Tx{tn) 



x(t n ) + (r - x{t n )) exp(-2TAi) 
where At = t n+1 - t n . Since x{t n ) = l^™! 1 / 2 and a:(i„ + i) = \ip n+1 \ 1/2 , it follows that 



(3.36) |V 



n+ i,_ ^ 1/2 lV-"l 



[IVI 2 + (r - |V'T) exp(-2rAi)] 1 /2 ' 
The phase 6 of ?/> is constant in time. If we multiply both sides by e l< ^, we obtain the 



expression (3.23) for the "semigroup" S. 



4. Evaluation. We now present the results of several experiments, where the 
algorithms described in the preceding section were applied to a benchmark problem. 

4.1. Benchmark Problem. The benchmark problem adopted for this investi- 
gation was the equilibration of a vortex configuration in a superconductor (Ginzburg- 
Landau parameter k = 16) embedded in a thin insulator (air), where the entire system 
was periodic in the direction of the free surfaces (y). 

The superconductor measured 128£ in the transverse (x) direction. The thickness 
of the insulating layer on cither side was taken to be 2£, so the total width of the 
system was 132£. The period in the y direction was taken to be 192£, so the entire 
configuration measured 132£ x 192£. 

The computational grid was uniform, with a mesh width h x = h y = i£. The 
periodic boundary conditions in the y direction were handled through ghost points, so 
the computational grid had 264 x 386 vertices. The index sets for the superconductor 
and blanket (see Eqs. (|2.16) and (|2.17)) were 



(4.1) 


Sc = {(i,j) :i = 5,.. 


. ,260,j = 


= 1,... ,386}, 


(4.2) 


Bl = {(t,j) :* = !,.. 


.,4,261,. 


..,264, j = l 



,386}. 

The applied field was uniform in y and equally strong on the left and right side of the 
system, 

(4.3) H L = H R = H = 0.5. 

(Units of H are H c ^/2, so H s=s 0.707 . . . H c ). As there is no transport current in the 
system, the solution of the TDGL equations tends to an equilibrium state. 

4.2. Benchmark Solution. First, preliminary runs were made to determine, for 
each algorithm, the optimal number of processors in a multiprocessing environment. 



Figure 4.1 shows the elapsed (wall clock) time for 50 time steps against the number of 
processors on the IBM SP2. Each algorithm showed a saturation around 16 processors, 
beyond which any improvement became marginal. All problems were subsequently 
run on 16 processors. 

Next, the fully explicit Alg orithm I was used to establish a benchmark equilibrium 



configuration. Equations (3.1)-(3.2) were integrated with a time step At = 0.0025 
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Fig. 4.1. Elapsed time for 50 time steps as a function of the number of processors. 



(units of £ 2 /D), the maximal value for which the algorithm remained stable. The 
evolution of the vortex configuration was followed by monitoring the number of vor- 
tices and their positions. Equilibrium was reached after 10,000,000 time steps, when 
the number of vortices remained constant and the vortex positions varied less than 
l.Ox 10~ 6 (units of £). The equil ibri um vortex configuration had 116 vortices arranged 
in a hexagonal pattern; see Fig. 4.2. The elapsed time for the entire computation was 
50.81 hours. The elapsed time per time step (0.018 seconds) is a measure for the 
computational cost of Algorithm I. 

4.3. Evaluation of Algorithms II— IV. Once the benchmark solution was in 
place, each of the remaining algorithms (II-IV) was evaluated for stability, accuracy, 
and computational cost. 

The stability limit was found by gradually increasing the time step and integrating 
until equilibrium. Above the stability limit, the algorithm failed because of arithmetic 
divergences. Equilibrium was defined by the same criteria as for the benchmark so- 
lution: no change in the number of vortices and a var iatio n in the vortex positions of 
less than 1.0 x 10 -6 . The results are given in Table 4T; At is the time step at the 
stability limit (units of £ 2 /D), N the number of time steps needed to reach equilib- 
rium, T the elapsed (wall clock) time (in hours) needed to compute the equilibrium 
configuration, and C the cost (in seconds per time step, C = 3600T/N). 

Because each algorithm defines its own path through phase space, one cannot ex- 
pect to find identical equilibrium configurations nor equilibrium configurations that 
are exactly the same as the benchmark. The equilibrium vortex configurations for the 
four algorithms were indeed different, albeit slightly. To measure the differences quan- 
titatively, we computed the following three parameters: (1) the number of vortices 
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Fig. 4.2. Equilibrium vortex configuration for the benchmark problem. 



in the superconducting region, (2) the mean bond length joining neighboring pairs of 
vortices, and (3) the mean bond angle subtended by neighboring bonds throughout the 
vortex lattice. In all cases, the number of vortices was the same (116); the mean bond 
length varied less than 1.0 x 1CP 3 £, and the mean bond angle varied by less than 
1.0 x 10~ 3 radians. Within these tolerances, the equilibrium vortex configurations 
were the same. 

Table 4.1 
Performance data for Algorithms I-IV. 





At 


N 


c 


T 


Algorithm 










I 


0.0025 


10,000,000 


0.018 


50.81 


II 


0.0500 


500,000 


0.103 


14.32 


III 


0.1000 


250,000 


0.232 


16.11 


IV 


0.1900 


131,580 


0.233 


8.41 



Finally, we evaluated the fully implicit Algorithm IV from the point of view of 
parallelism. From the benchmark problem we derived two more problems by twice 
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doubling the size of the system in each direction, while keeping the mesh width the 
same (tj£). The resulting computational grid had 528 x 772 vertices for the intermedi- 
ate problem and 1056 x 1544 vertices for the largest problem. Speedup was defined as 
the ratio of the wall clock time (exclusive of I/O) to reach equilibrium on p processors 
divided by the time to reach equilibrium on a single processor for the smallest and 
intermediate problem, or twice the time to reach equilibrium on two processors for the 
largest problem. (Th e largest problem did not fit on a single processor.) The results 
are given in Fig. 4.3. The curve for the benchmark problem was obtained as an aver- 









-V ' 


s® ~° - 




^—1056x1544 
-a- 528x772 
- .-264x386 


















f 


_y / 




_--«-" 








r 


s* 












/ 

















10 15 20 25 30 35 40 

Number of Processors 



Fig. 4.3. Computational cell with evaluation points for i/>, A x , and A y . 

age over many runs; the data for the intermediate and largest problem were obtained 
from single runs, hence they are less smooth. The speedup is clearly linear when the 
number of processors is small; it becomes sublinear at about 12 processors for the 
smallest problem, 14 processors for the intermediate problem, and 18 processors for 
the largest problem. 

5. Conclusions. The results of the investigation lead to the following conclu- 
sions. 

(1) One can increase the time step At nearly 80-fold, without losing stability, by 
going from the fully explicit Algorithm I to the fully implicit Algorithm IV. 

(2) As one goes to the fully implicit Algorithm IV, the complexity of the matrix 
calculations and, hence, the cost C of a single time step increase. 

(3) The increase in the cost C per time step is more than offset by the increase 
in the size of the time step At. In fact, the wall clock time needed to compute the 
same equilibrium state with the fully implicit Algorithm IV is one-sixth of the wall 
clock time for the fully explicit Algorithm I. 

(4) The (physical) time to reach equilibrium — that is, NAt, the number of time 
steps needed to reach equilibrium times the step size — is (approximately) the same 
for all algorithms, namely, 25,000 (units of £ 2 /D). 

(5) The fully implicit Algorithm IV displays linear speedup in a multiprocess- 
ing environment. The speedup curves show sublinear behavior when the number of 
processors is large. 
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